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Abstract 

In this paper we propose a highly accurate approximate performance analysis of a het¬ 
erogeneous server system with a processor sharing (PS) service discipline and a general 
job-size distribution under a generalized join the shortest queue (GJSQ) routing protocol. 
The GJSQ routing protocol is a natural extension of the well-known join the shortest queue 
(JSQ) routing policy that takes into account the non-identical service rates in addition to 
the number of jobs at each server. The performance metrics that are of interest here are 
the equilibrium distribution and the mean and standard deviation of the number of jobs at 
each server. We show that the latter metrics are near-insensitive to the job-size distribution 
using simulation experiments. By applying a single queue approximation (SQA) we model 
each server as a single server queue with a state-dependent arrival process, independent 
of other servers in the system, and derive the distribution of the number of jobs at the 
server. These state-dependent arrival rates are intended to capture the inherent correlation 
between servers in the original system and behave in a rather atypical way. 


1 Introduction 

1.1 Motivation 

This work is motivated by web server farms. Server farms have gained popularity for pro¬ 
viding scalable and reliable computing and web services. Most commonly the objective in 
analyzing such a system lies in the determination of an optimal or near-optimal load balancing 
routing protocol so as to maximize the performance of the system, see, e.g., HlEli, where 
the performance of interest is usually the mean response time for an arbitrary job. In this 
paper the objective is to report some interesting properties of the arrival flow to each server 
and suggest an approximation approach for the GJSQ routing protocol. We consider farms 
with heterogeneous servers, which is motivated by the different hardware and the wide variety 
of computing capacities regarding processing power and memory access performance seen in 
practice in server farms m- We assume that service requests arrive to the system according 
to a Poisson process. Upon arrival, a front-end dispatcher routes the request to one of the 
servers. After the request has been routed to the server, we assume that it cannot balk or 
jokey. All requests routed to a server are sharing the provided service (think of bandwidth, 
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CPU, or RAM). We assume a PS service discipline at each server since it closely approximates 
the scheduling policies [zmn] employed by most commodity operating systems (e.g., Linux 
CPU time-sharing) and is a popular policy in computing centers (e.g., Cisco Local Director, 
IBM Network Dispatcher and Microsoft Sharepoint, see [3] for a survey). 

In [5] the authors consider a server farm consisting of homogeneous servers, where upon 
arrival jobs are routed according to the JSQ routing protocol. This protocol in case of homo¬ 
geneous servers, due to the PS service discipline, is performing near-optimal in terms of the 
mean response time. However, as indicated by Whitt in |15j . the JSQ policy is far from optimal 
in case of heterogeneous servers. In [2] the authors comment on the performance of various 
systems under different routing protocols and conclude that the shortest expected delay (SED) 
routing protocol is near-optimal in terms of mean response time. The SED policy is a policy 
that routes jobs upon arrival to the queue promising the minimum expected delay (which also 
includes the processing time). In case of exponential job-size distributions, the GJSQ and SED 
routing protocols are identical and in case of homogeneous servers GJSQ and JSQ are the 
same. However, in case of general job-size distributions and heterogeneous servers we assume 
that the only available information are the service rates and the number of jobs at each server, 
i.e. we do not keep track of residual processing times. 

Due to the complexity and the various challenges that the model at hand presents, we 
restrict our analysis to the case of two heterogeneous servers with a general job-size distri¬ 
bution under the GJSQ routing protocol. Erom here onwards we refer to this model as the 
M/G(l, s)/2/GJSQ/PS system, abbreviated as the GJSQ model, where G is the job-size dis¬ 
tribution and 1 and s are the service rates at servers 1 and 2, respectively. The approach 
described in this paper can be seen as a first stepping stone towards the analysis of heteroge¬ 
neous server farms with PS servers; a very broad area, full of interesting problems. Moreover, 
the ideas presented here extend the work of Gupta et al. [5] on the analysis of the JSQ routing 
for homogeneous web server farms. 

1.2 Related work 

To the best of our knowledge there is no previous mathematical analysis of the GJSQ system. 
In [13], Selen et al. derive the joint equilibrium distribution of the number of jobs at each server 
in the M/M(l, s)/2/SED/FGFS model. They prove that this distribution can be expressed 
as an infinite series of product forms using the compensation approach. The benefit of that 
approach is that it produces, by truncating the series expression, numerical results with an a 
priori set accuracy level. Unfortunately, the compensation approach is not appropriate (in its 
current setting) for multiple servers, nor for general job-size distributions. Before |l3|, very little 
was known regarding the mathematical analysis of the SED policy. In m , the authors suggest 
two models that act as upper and lower bounds to the SED system. However, they do not 
provide closed form expressions for the equilibrium distribution of these two bounding models, 
but only an algorithmic approach based on matrix analytic methods. Furthermore, in mi, 
the authors show that the SED routing policy is asymptotically optimal in terms of the mean 
response time and results in complete resource pooling in the heavy traffic limit. This heavy 
traffic limit result may be used in a similar manner as in |12j . However, after a few numerical 
experiments, we concluded that this approximation in our case results in poor estimates and 
for this reason we did not proceed in this direction. On the contrary, the approach developed 
by Gupta et al. |5| on approximating the distribution of the number of jobs at each server, 
as we show in this paper, is appropriate for the GJSQ setting with heterogeneous servers. 
More concretely, in |5], the authors develop the SQA method that accurately determines the 
distribution of the number of jobs at each server by modeling each queue as an Mn/M/1/ PS 
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system with state-dependent arrival rates. These state-dependent arrival rates are referred 
to as the conditional arrival rates and are constructed in such a way that they capture the 
inherent correlated behavior of the complete server farm. 


1.3 Contributions 

We believe that we provide the first approximate analysis of the equilibrium distribution and 
moments of the number of jobs at each server in the GJSQ system (and by Little’s law also the 
mean response time for an arbitrary job). Moreover, the approximation is highly accurate: we 
encounter a maximum relative difference between the approximation and simulations of 2.2%. 
In deriving these approximations, we provide three key contributions: 


1. The mean and standard deviation of the number of jobs at each server and the conditional 
arrival rates are near-insensitive to the job-size distribution. This allows us to study the 
more tractable model with an exponential job-size distribution. 

2. In case of an exponential job-size distribution, the SQA method yields the same equilib¬ 
rium distribution for the number of jobs at each server as in the original GJSQ model. 

3. For the application of the SQA method we present an approach for the derivation of the 
conditional arrival rates. In particular, we show that the conditional arrival rates, say 
Aj(re), i = 1, 2, n G No, to server 1 satisfy 

Ai(n) —)• as n —>■ oo, (1-1) 


where p is the load on the system, see Section and the conditional rates to server 2 for 
large n oscillate between s different points. Note that the former result is similar to the 
result obtained in [5] for the case s = 1, however the latter result is very atypical and is 


discussed in greater detail in Section 2.3 


1.4 Outline 

The rest of the paper is organized as follows. In Sectionwe give a detailed model description 
and formally define and investigate the time-average and conditional arrival rates. Section is 
devoted to showing that the performance metrics of interest are near-insensitive to the job-size 
distribution. We describe the SQA and determine the conditional arrival rates in Section 
The approximations are evaluated in Section In Section we present some conclusions. 


2 Model description 


2.1 Heterogeneous servers 

We consider a system of two heterogeneous servers and a single dispatcher. The servers employ 
a PS service discipline and can have different service rates, i.e. server I has service rate 1 and 
server 2 has service rate s. Jobs arrive to the dispatcher according to a Poisson process with 
rate A and are routed immediately to one of the servers. Jobs cannot switch servers after being 
routed. We detail the routing policy in Section 2.2 The size of a job is drawn from a general 


distribution G. Without loss of generality we assume that the mean job size is 1. Note that, 
for example, the (residual) processing time of a (residual) size G job that runs on server 2 that 
is currently serving q 2 jobs is given by Gq 2 /s. 
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Figure 1: Simulated long-term fraction of jobs routed to server i as a function of the load 
p, where s = 4 and the job-size distribution is exponential. Dashed lines represent expected 
behavior. 


In what follows we assume that s is a positive integer number. In the general case s G M+ 
we can bound the corresponding system by two systems with service rates given by the closest 
two integers to s. 

2.2 Routing policy 

The routing policy employed by the dispatcher is a state-aware policy, i.e. the dispatcher is 
aware of the number of jobs at each server just before an arrival instant, qi and q 2 , and the 
service rates. The GJSQ routing policy routes an arriving job to the server with the smallest 
index (qi -|- l)/si, where Si is the service rate at server i. In case of a tie, the job is randomly 
routed to one of the servers. These indexes may be interpreted as an estimate of the expected 
processing time for the arriving job, made by the dispatcher who is unaware of the job-size 
distribution and the remaining processing times of the jobs currently in service, and furthermore 
ignores future arrivals. 

Under this routing policy, we define the load on this system as 

p:=X/{l + s). (2.1) 

Throughout the rest of this paper we assume that p < 1. 

Although not necessarily optimal, GJSQ routing outperforms JSQ routing when servers are 
non-identical. GJSQ routing attempts to balance the load on the servers by taking into account 
the different service rates in addition to the information on the current number of jobs at each 
server. In Figure we show that the long-term fraction of jobs routed to the two servers is a 
function of the load p. In light traffic GJSQ assigns all jobs to the fast server and in heavy 
traffic the load is divided proportionally according to the service rates. This is in contrast with 
JSQ routing, which assigns a long-term fraction of the jobs to server 1 that decreases from 1/2 
to 1/(1 -|- s) for increasing load p (verified through simulation). 

2.3 Arrival rates 

We briefly introduce two important concepts related to the (time-average) arrival rates to each 
server. These concepts will be used throughout the paper. 

Definition 2.1. In the GJSQ model, the time-average arrival rate to server i is defined as 

A := lim (2.2) 

t^oo t 
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Figure 2: The conditional arrival rates to server 2 oscillate between s points. 


where Ai(t) is the number of arrivals at server i during the time interval [0,t]- 


Definition 2.2. In the GJSQ model, the conditional arrival rate to server i, given that server 
i has n jobs, is defined as 


Aj(n) := lim 

t^oo 


Ti,n{t) ’ 


(2.3) 


where Ai^n{t) is the number of arrivals at server i during the time interval [0, t] that see n jobs 
at server i on arrival (excluding themselves), and is the total time spent by server i with 

n jobs during the time interval [0,t]. 

The two definitions above are related. Assuming it exists, let 7rj(n) be the equilibrium 
probability that there are n jobs at server z, then A* = -^*(^)'^i(’^)- 

Figure [^depicts the conditional arrival rates to server 2 for varying s. Intuitively it makes 
sense that if a server has many jobs, the other server will probably have few jobs and thus it 
is less likely that the dispatcher routes the job to that server. However, what we see here is 
a peculiar repeating pattern that has s different points and does not align with this intuition. 
We see that if server 2 has a multiple of s jobs (or one less), fewer jobs are routed to server 2. 
This pattern is difficult to explain, but it is definitely related to the probability that server 1 
has a lower index than server 2, given that server 2 currently has n jobs. We expect and indeed 
verify that this probability also follows such a repeating pattern. Additionally, states in server 
2 are somewhat similar if they differ by a multiple of s jobs, which can be derived from the 
equilibrium distribution in |14) . 


3 Near-insensitivity 

In [2j the authors establish a near-sensitivity property in the setting of a homogeneous server 
farm with JSQ routing. In particular, the first and second moment of the number of jobs at 
server i, Qi, and the conditional arrival rates are near-insensitive to the job-size distribution. 
The near-insensitivity of these two metrics seems related to the insensitivity of the equilibrium 
distribution to the job-size distribution in PS servers, see, e.g., [5l Theorem 4.2]; and the fact 
that the routing policy only uses the number of jobs at each server when making a decision, 
as opposed to, e.g., using residual processing times. The GJSQ routing decisions are based on 
the dynamically changing number of jobs at each server as well as the service rates. Indeed, 
one expects the near-insensitivity properties to extend also to the case of heterogeneous servers 
and GJSQ routing. Establishing this near-insensitivity property is important, since it allows 
us to limit our attention to the more tractable GJSQ system with an exponential job-size 
distribution. 
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Name 

Distribution 

Support 

Variance 

uni 

Uniform 

[0,2] 

1/3 

exp 

Exponential 

[0,oo) 

1 

weib 

Weibull 

(0,oo) 

5 

logn 

Log-normal 

(0,oo) 

10 


Table 1: Job-size distributions used in simulations. 


Job-size distribution SQA 


s p 

Metric 

uni 

exp 

weib 

logn 

Value 

Dili. 

2 0.7 

E[Qi] 

0.9139 

(0.0030) 

0.9232 

(0.0030) 

0.9223 

(0.0029) 

0.9361 

(0.0038) 

0.9077 

1.7% 


cr(Qi) 

1.0404 

(0.0049) 

1.0505 

(0.0050) 

1.0560 

(0.0056) 

1.0704 

(0.0067) 

1.0462 

0.4% 


E[02] 

2.0111 

(0.0059) 

2.0289 

(0.0061) 

2.0222 

(0.0057) 

2.0519 

(0.0074) 

2.0329 

-0.2% 


(j{Q2) 

2.0302 

(0.0099) 

2.0465 

(0.0106) 

2.0506 

(0.0114) 

2.0813 

(0.0141) 

2.0484 

-0.1% 

0.9 

E[Qi] 

3.2244 

(0.0336) 

3.2797 

(0.0336) 

3.2316 

(0.0298) 

3.2396 

(0.0266) 

3.2188 

1.9% 


cr(Qi) 

3.2208 

(0.0590) 

3.2716 

(0.0723) 

3.2186 

(0.0575) 

3.2002 

(0.0505) 

3.2161 

1.7% 


E[Q2] 

6.6841 

(0.0676) 

6.7915 

(0.0674) 

6.6834 

(0.0587) 

6.6988 

(0.0524) 

6.6424 

2.2% 


cr(Q2) 

6.4289 

(0.1185) 

6.5288 

(0.1453) 

6.4186 

(0.1141) 

6.3828 

(0.1016) 

6.4091 

1.9% 

4 0.7 

E[Qi] 

0.4688 

(0.0016) 

0.4747 

(0.0017) 

0.4705 

(0.0017) 

0.4667 

(0.0022) 

0.4741 

0.1% 


cr(Qi) 

0.6685 

(0.0024) 

0.6730 

(0.0026) 

0.6700 

(0.0029) 

0.6652 

(0.0031) 

0.6655 

1.1% 


E[Q2] 

2.5386 

(0.0063) 

2.5507 

(0.0069) 

2.5177 

(0.0070) 

2.4997 

(0.0067) 

2.5866 

-1.4% 


cr(Q2) 

2.5082 

(0.0102) 

2.5179 

(0.0115) 

2.4936 

(0.0133) 

2.4744 

(0.0122) 

2.5457 

-1.1% 

0.9 

E[Qi| 

1.8662 

(0.0191) 

1.8793 

(0.0145) 

1.8830 

(0.0196) 

1.9400 

(0.0223) 

1.8813 

-0.1% 


cr(Qi) 

1.9404 

(0.0338) 

1.9539 

(0.0314) 

1.9801 

(0.0444) 

2.0394 

(0.0408) 

1.9566 

-0.1% 


E[02] 

8.2405 

(0.0769) 

8.2773 

(0.0597) 

8.2631 

(0.0783) 

8.4863 

(0.0861) 

8.3642 

-1.0% 


cr(Q2) 

7.6982 

(0.1374) 

7.7507 

(0.1264) 

7.8485 

(0.1815) 

8.0912 

(0.1652) 

7.7692 

-0.2% 


Table 2: Simulated mean and standard deviation of Qi, for the GJSQ system with various 
s, p and job-size distribntions. Sample standard deviation is shown in parentheses. Last two 
columns show the value obtained by the SQA and the relative difference with respect to the 
exponential case. 


3.1 Simulation settings 

To support our claims, we simulate the GJSQ model. A simulation consists of 2 • 10® job 
departures and each simulation is repeated 50 times. Statistics are only computed for departed 
jobs, i.e. data of jobs that are still in service at the end of the simulation are discarded. In 
Table we list the four job-size distributions considered in this paper. 

3.2 Near-insensitivity results 

In Table we show simulated statistics on the mean and standard deviation cr(-) of Qi for 
the GJSQ model with various job-size distributions. For the settings considered in Table 
the mean nnmber of jobs at server i deviates by no more than 3.6% from the exponential 
case, while the standard deviation deviates by at most 4.4%. The largest deviations from the 
exponential case occur for the log-normal job-size distribntion. This is as expected, since this 
job-size distribntion has a variance that is 10 times higher than the variance of the exponential 
job-size distribntion. Although the results are not as strong as those shown in O Figure 3], we 
conclude that the more volatile environment of heterogeneous servers and GJSQ routing also 
has the near-insensitivity property for E[(5i] and (j{(Qi). Moreover, the performance in terms 
of the mean response time is also near-insensitive to the job-size distributions by Little’s law. 
Goncerning the conditional arrival rates, we see in Fignre that the simnlated values for 
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n n 

(a) p = 0.7, s = 4, A = 3.5 (b) p = 0.9, s = 4, A = 4.5 


Figure 3: Simulated conditional arrival rates in the GJSQ system with various job-size dis¬ 
tributions. The dotted curves represent values determined by the algorithm in [T3] for the 
exponential job-size distribution. 


the job-size distributions of Tablematch the results of the algorithm for the exponential case 
M- Simulated values for states where the sample standard deviation is not too high differ by 
at most 5% from the results for the exponential case. So, also the conditional arrival rates are 
near-insensitive to the job-size distribution. 


4 Single queue approximation 


We have established near-insensitivity of E[Qj], cr{Qi) and the conditional arrival rates to the 
job-size distributions. Thus, we may limit our attention to systems with an exponential job- 
size distribution. In this section we derive an approximation for the distribution of the number 
of jobs at each server using the SQA, which models server i as an Mn/Mi/1/PS queue with 
state-dependent arrival rates Aj(n), see also [5l Section 3]. SQA is exact when the job-size 
distribution is exponential and the routing belongs to a specific class of routing policies; the 
following theorem is a version of [5l Theorem 3.1] that is applicable to the GJSQ model. 

Definition 4.1. A stationary state-aware routing policy is a time-stationary routing policy 
that only uses information about the number of jobs at the servers and the service rates at the 
instant of an arrival. The decisions may be made probabilistically, possibly biased in favor of 
certain servers. 


Theorem 4.2. Consider the MlM{\,s)l2lTZlS queueing model, where TZ is any station¬ 
ary state-aware routing policy, e.g. GJSQ, and S is any stationary, size-independent, work- 
conserving service discipline, e.g. PS. Consider server i in this model. Then SQA with the 
exact conditional arrival rates Ai(-) yields the same equilibrium distribution for the number of 
jobs at each server as in the original model. 


It remains to specify the conditional arrival rates Aj(n) for both servers. We combine exact 
limiting results for n > Ni and approximation results for n < W, where Ni = 3 and N 2 = 2s. 
These choices for W result in accurate approximations. 

We note that Theorem 4.2 implies that in order to determine the conditional arrival rates, 
we may assume a FGFS service discipline. 

In Figure [^we have seen that the conditional arrival rates Aj(n) exhibit a repeating pattern 
from some n and onwards. We rigorously characterize this limiting repeating pattern in the 
next theorem. 
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Theorem 4.3. For the M/M(l, s)/2/GJSQ/PS queueing model with s G N, 


^hm Ai(n) = 

n^oo 


A 2 ™(r) := lim A 2 (sn + r) = 


* A{r) ’ 
,,1+s ^(0) 
pP A(s-1)> 


r = 0,1,... ,s - 2, 
r = s — 1, 


(4.1) 

(4.2) 


where 

A{r) = + /i(^) + (4.3) 

and the variables /3i,/32, • • • ,/3s+i, ^ 1 ,^ 2 , ■ ■ ■ ,Ps, o.nd the functions h{-), i+(-)j *-(•) defined 
in Appendix P] 

Proof. See Appendix [B| □ 

For the rates Ai(n), n < 3 and A 2 (n), n < 2s we provide approximations that are functions 
of s and p. For server 1 we use a multiple linear regression model to fit an approximate function 
for the conditional arrival rates on data obtained from the algorithm in [I4j for s = 1,2,3,4 
and p from 0.3 to 0.99. Obviously, one can also use conditional arrival rates obtained by 
simulation for these fitting purposes. We carefully select a set of 5 independent variables for 
each conditional arrival rate. This leads to the following approximate conditional arrival rates 
for server 1: 


MO) 

pl+s 

Ml) 

pl+s 

Ai(2) 

pl+s 


sp 

^ 7 

1 4 

[0. 

669 

-1.90 

1.23 

1.86 -0.192]'^, 


sp^ 

1 1 
p 

± filp 

sp ^ 

[-0 

.00856 

1.37 

-0.0578 0.123 

-0.254]^, 

1 + 

sp i 

p 

p 

1 

7 ^ 

pi/s 

— [- 
100 L 

-0.131 

-0.820 -6.48 

10.4 0.893]^ 


(4.4) 

(4.5) 

(4.6) 


with the transpose of a vector x. For s = 1, one should consider Ai(-) = A 2 (-) and use the 
approximations presented in (4.4)-(4.6). 

For server 2, let us note that A 2 (n) = A, n = 0,..., s — 2 due to the GJSQ routing. Using a 
multiple regression model in this case is more difficult, since the number of states for which we 
need to obtain a fit increases with s. To circumvent a possibly complex fitting procedure, we 
establish a relation between the conditional arrival rates for the states n = s — l,s,...,2s — 1 
and the limiting conditional arrival rates determined in Theorem 2. Namely, 


A2(n) ~ (1 + 


2s - 


l) 


1 

2 n-{s-l) 



s), 


(4.7) 


behave in various 

limiting regimes as expected: 

Proposition 4.4. 


where for convenience A 2 ™(—1) = A 2 ™(s —1). The approximations (4.4)-(4.7 


1. For s -+ 00 , we have that Ai(n) f 0 and \ 2 {n) = A for all n G Nq. No job will join server 
1, since the processing times in server 2 are instantaneous. 


2. In the light-traffic regime, i.e. p f 0, we find that Ai(n) | 0, n G Nq and A 2 (n) i 0, n > 
s — 1. 
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n n 

(a) s = 3, server 1 (b) s = 3, server 2 



(c) s = 4, server 1 


(d) s = 4, server 2 


Figure 4: Comparison of conditional arrival rates determined by the algorithm in [T3] (lines) 
and our approximations (marks) for p = 0.4 (—,A), p = 0.7 (--,o), and p = 0.9 ( ,□). 


3. In the heavy-traffic regime, i.e. p \ 1, we establish that Ai/A = 1/(1 + s) and A 2 /A = 
s/(l + s) which is consistent with the findings in Figure^ 

Proof. 1. Follows straightforwardly by taking the limit s —)• 00 in ( |4.4[ )-(4.6 ) while taking into 
account that p = A/(l + s). Furthermore, observe that X 2 {n) = A, n = 0,..., s — 2, so that 
lims^oo A 2 (n) = A, n G Nq. 

2. See Appendix [C| 

3. From the approximate conditional arrival rates Ai(-) one can derive (approximate) equi¬ 
librium probabilities 7ri(-). Then, Ai = ~ -|- 1) = 1 — vri(O) 

by exploiting the balance equations. For p f 1 it can be verified that 7ri(0) f 0, so that 
limp'i'i Ai/A = 1/(1 -|- s). The result for server 2 follows analogously. □ 


5 Evaluating the approximation 

We are now in a position to evaluate the proposed approximations. First, we show that 
the approximations for the conditional arrival rates follow closely the exact values, which were 
determined using the algorithm |14] . Second, we establish that the mean and standard deviation 
of the number of jobs at each server is also well approximated. 

Figure compares the conditional arrival rates obtained from the algorithm in [TTj and 
the approximations derived in the previous section. For the cases considered in the hgure, the 
maximum relative difference of the approximation with respect to the values determined by the 
algorithm is 1.5% for Ai(-) and 4.1% for A 2 (-)- Since both methods consider exponential job-size 
distributions, the difference is due to the fitting error introduced in the approximations of the 
conditional arrival rates in Section]^ and the truncation error in the algorithm in m. However, 
since the truncation error has been chosen to be of the order 10“®, it has little influence. 

In the two rightmost columns of Table we provide the mean and standard deviation of the 
number of jobs at both servers determined using the SQA. We report highly accurate approxi¬ 
mations that deviate less than 2.2% from the case with an exponential job-size distribution for 
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Figure 5: Simulated conditional arrival rates for a system with three servers with service rates 
1 ( ), 2 and 5 (—), with p = 0.7. 


the listed values of s and p. Although our approximations are not aimed at the case s = 1, we 
report accurate approximations also in this setting with maximum relative differences of the 
same order as in 0 Section 6.1]. 


6 Conclusion 


In this paper, we provide an approximate performance analysis of a queueing system consisting 
of two heterogeneous PS servers with service rates 1 and s G N, respectively, a general job- 
size distribution and GJSQ routing. More concretely, we derived the approximate equilibrium 
distribution of the number of jobs at each server using the SQA method. In order to apply 
SQA, we established that the GJSQ system is near-insensitive to the job-size distribution 
and thus we approximated the system at hand with exponentially distributed job-sizes. We 
then approximated the conditional arrival rates for the exponential case, by combining exact 
limiting results for large number of jobs and approximation results, which were obtained using 
a multiple linear regression model, for small number of jobs. Ultimately, the aforementioned 
approach resulted in approximations that are highly accurate; we reported a maximum relative 
difference with respect to exact or simulation results of 4.1% for the conditional arrival rates 
and 2.2% for the mean and standard deviation of the number of jobs at each server. 

In this paper we set the groundwork for the analysis of server farms with heterogeneous 
servers under the GJSQ routing policy by analyzing the case of two servers. Of course, server 
farms consist of multiple servers so it is in our future plans to extend the analysis presented 
in this paper to more than two servers. The most difficult aspect of this task would be the 
derivation of the conditional arrival rates, which possibly has to rely on simulation data, since 
the approach in m is in its current setting restricted to two servers. In Figure]^ we present an 
example of the simulated conditional arrival rates in case of three servers with service rates I, 
2 and 5. Note that the structure of the conditional arrival rates is as expected, i.e. the number 
of points in the repeating pattern is directly related to the rate of service, but the exact values 


of these points differ from the values obtained by formulas (4.1) and (4.2). 
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A Definition of the variables and functions in Theorem 4.3 


The definitions can be found in HI Lemma 5.11], but we summarize them here. We denote 


a = 


The functions i+(-) and i_(-) are vectors of size s and have their first element equal to 1, 
i.e. i+(a,/3,0) = i_(a,/3,0) = 1. Furthermore, the vectors satisfy 


z+(a,/3,r) /q;/ 3(1 + s)(/5+ 1) - /3^(1 + s)p - 


f+(a,/5,0) 


aj3s 


) , r = 0,1,..., 


s — 1 


(A.l) 


and 


i_ia,l3,r) ^ - T(a,/3, V;+(/3))V;_(/3)^ 

z_(a, /3,0) 


, r = 0,l,...,s-l (A.2) 


with 


and 


y±iP) = 


T(a,/3, V'-C/^)) - T(q;,/3,V'+(/3)) 

{1 + s){p + 1) - (5 ± - {1 + s){p + l)y - 4s(l + s)p 


2s 


(A.3) 


^{a, = P - {1 + s){p + 1) + sy +-{I + s)p'tp'’ ^ = {l + s)p'4> V-V^®-l). (A.4) 

a a 

The variables j3i, /32, ■ ■ ■, Ps are the s roots that satisfy |/3j| < ja], i = 1,2, ...,s of the 
equation 

(a/3(l + s)(/j + 1) - ,5^(1 + s)p - a'^y - (3{a(3sY = 0, (A.5) 

and f3s+i with |/3s+i| < jaj is the single root of 


+ /3^((i + s)py - aysyy+iyy + y-iyy) = o. 

The vector h = (/i(0), h{l ),..., h{s — 1)) is given by 


h = a(^(l + ^^+al) ^'^pii+{a, Yi), 


(A.6) 


(A.7) 


i=l 


where the coefficients pi,p 2 , ■ ■ ■ ,^3 satisfy 

S 

rji {jii ((1 + s)pl + 

i=l 

+ (—(1 + s )(/9 + 1 )/+ + (1 + —(1 + + a /]) ^ i _|.( Q :,/ 3 j ) 

= -/3^+i((l + + al)\_{a, I3s+i), (A.8) 

where I is the sx s identity matrix, is the sx s binary matrix with element (x, y) equal 

to one and zeros elsewhere, and L is an s x s subdiagonal matrix with elements (x, x — 1), x = 
1, 2,..., s — 1 equal to one and zeros elsewhere. For consistency with indexing of all vectors, 
the indexing of a matrix starts at 0. 
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4.3 


B Proof of Theorem 

The proof is based on the exact results of the related M/M(l, s)/2/SED/FCFS system, with 
s G N, presented in jlTj. Although we obtain similar results for the limiting conditional arrival 
rates for server 1 as in [5], we use here a completely different approach in deriving the limits. 

In [2], the state space {{qi,q 2 ) \ {qi,q 2 ) £ Nq} Markov process is transformed to 

the state space {(m,n, r) | m G No, n G Z, r = 0,1,..., s — 1} where m = min(gi,[^J), 
n = — qi and r = mod(g 2 ,s). Let us denote the equilibrium probabilities for the three- 

dimensional state space as p{m,n,r). The equilibrium probability p{m,n,r) has a series ex¬ 
pression, i.e. p{m,n,r) = namely, for m > 0, n > 1, 

OO (•5 + 1)* s 

p{m,n,r) = CY^ ^ 

1=0 i=l j=l 


For m > 0, 


OO (-5+1)^ 

p{m,0,r) = CY^ af^h^ir). 


1=0 i=l 


(B.2) 


For m > 0, n < —1, 


OO (s+1)' 

p{m,n,r) = CY («;,*, A,i(^+i), 0 

1=0 1=1 


OO (s+1)' 

/=0 i=l 


(B.3) 


For the exact interpretation of each variable we refer the reader to [13]. In m the authors 
establish the following properties: 


1. There exists a positive integer N such that the series in (B.l) and (B.3) converge abso¬ 
lutely for all m > 0, |n| > 1 with m -|- |n| > N and the series (B.2) converges absolutely 
for all m > N. 


2. For m -|- |n| > N, we have \x{l, m, n,r)\ < u{l) and '*^(0 < oo. 

3. The series n, r), r = 0,1,..., s — 1 converges absolutely. 

4. \ai^i\ > \lii,d{i)+j\ and |A,j| > |az+i,il with ao,i = < 1. 

In this proof we make use of the dominated convergence theorem for complex-valued functions. 


B.l Server 1 

The limiting conditional arrival rate to server 1 can be determined from 


lim Ai(n) = lim 

n—)-oo n^oo 


7ri(n + l)/ag|^ 
7ri(n)/a[fi 


ao,i- 


(B.l) 
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The marginal distribution for server 1 is given by, where m = [yj and r = mod(g 2 ; •s)) 


oo s—1 

7ri(n) = EE p(min(n, m),m — n, r) 

m=0 r=0 

n—1 s—1 s—1 oo s—1 

= EE p{m, m — n,r) + 0, r) + EE p{n, m, r). 

m=0 r=0 r=0 m=l r=0 


Furthermore, 


lim 

n^oo 


7ri(n) 

“o,i 


p{m, m — n,r) 


n—1 5—1 

: lim yy 

“ 0,1 

+ v lim + y; V lim Kwf, 

n-s>oo an 1 n-l-oo Qin i 

r=0 0,1 m=l r=0 0,1 


(B.5) 


(B.6) 


where the interchange of the limit and the series for the third term on the right-hand side of 


(B.6) is allowed by the dominated convergence theorem, because one can bound p(n, m, r) from 


above by p{0,m,r) and °° since it is a subseries of X]m+|n|>Af 

which converges absolutely by property 3. Furthermore, we know that limm-i-oo pirn, n, r) = 
limrn^oo which is equal to x{l, m, n, r) by the dominated con¬ 

vergence theorem for complex-valued functions in combination with property 2. This allows us 


to compute the second and third term on the right-hand side of (B.6). The first term on the 


right-hand side of (B.6) can be determined as follows 


n—1 s—1 , , 

,. p(m,m — n,rj 

hm > > —---- 

n—>-oo ^ 


m=0 r=0 


“0,1 


(.+1)' 


00,1 


n /a \n 

P^,^(3 + l) ^ 
OO 


Ol,i 


r=0 


E E "WO+D- ^_ 

/=0 i=l A,i(s+i) 

oo (s-l-1)' + / o^l + l.ga + l) 

+ ^ ^^+hh^+i) ^ . oi+i^(3+? --yi-(«i+i,i(,+i),A,*(s+i),r)) 

A,i{s+i) 


1=0 i=l 
:-l 


CJ- J. 

= C ^ y i_(«o,i,/3o,^+i,r). 

/9o,s+l r=0 


(B.7) 


Interchange of the limit and series is again allowed here since one can bound the absolute value 
of the summands from above by v{l) for sufficiently large n and ^ finally 

establish that 


lim llM = c 


n^oo a, 


0,1 


'>10,s+l 
ao 

/5o,s+i 
5—1 


5—1 




r=0 


r) 

5-1 


+ y ho,i{r) + y y *+(ao,i, /3o,j, ' 

r=0 7=1 


Thus, by (B.4) and (B.8), hm„_>.oo Ai(n) = ao,i = P 


— 


(B.8) 
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B.2 Server 2 

The limiting conditional arrival rate to server 2 can be determined from 

r = 0,l,...,s-2, 


lim X 2 {sn + r) = 

n—^oo 


^7r2(sn+r+l)/Q:g'j 
-■■n—^oo7r2 (sn+ri/cfp ’ 
^7r2(sn+r+l)/oQ’|^ 
7r2(sn+r)/a^^^ 


lim^^oo S «0,1 ^ r = s-l. 


(B.9) 


The marginal distribution for server 2 is given by, for r = 0,l,...,s — 1, 

/ N / , sn + r .. . sn + r . .. 

7r2(sn + r)= } pfmmfgi, I-I), I-I - gi, modfsn + r, sj) 

s s 

91=0 

OO 

= p(min(gi,n),n - qi,r) 

91=0 

n—1 OO 

= p{qi,n - qi,r) + p{n,0,r) + p{n,-qi,r). 

91=0 qi=l 


(B.IO) 


For 7T2{sn + r + l), r = s — Iwe should replace n by n + 1 and r by 0 in (B.IO). Furthermore, 
for r = 0,1,..., s — 1, 


lim = lim V + lim + lim V (B.n) 

77,—Ino ^. n^cta . t 7.—Ino -I—^ 


n^oo CK, 


0,1 


91=0 


a 


0,1 


n—loo Q; 


0,1 


91=1 


a, 


0,1 


Using identical arguments as for the limiting conditional arrival rate for server 1, we establish 

-K 2 {sn + r) 


lim 

n—>-oo a 


= A{r), r = 0,l,...,s-l. 


(B.12) 


0,1 


where A{r) is given in (4.3). Finally, combining (B.9) and (B.12) proves (4.2). 


C Proof of Proposition 4.4, point 2 


By letting /? 0 in ( |4.4[ )-( [4^ and Ai(n) « n > 3 we immediately find that Ai(n) 0, n G 

No. 

We note that in (4.7) the factors on the right-hand side in front of A 2 ™(u—s) go to a constant 
for p I 0. So, what remains is that we establish that limp^o A 2 ™(i’) = 0, r = 0,1,..., s — 1. 
This part of the proof relies heavily on the asymptotic results of [13] . We denote a = and 
investigate for r = 0,1,..., s — 1, 


A{a,r) ^ Pi/a /Pisr/s , lu 


a 


r/s 


i=l 


1 - Pi/a 


h{r) 


a 


Q,(r+l)/s 


-|- a 


l-r/s Ps+l/ct ^ 


1 - Ps+1 


i.{a,Ps+i,r), (C.l) 


where we used that i+(Q:, Pi, r) = UiP[^^ with Ui the f-th unit root of = 1, which is established 
in m Lemma 5.6]. Now, 


A{a,r) 

iim-;— = c(r), 


(C.2) 


alo a’’A 

where c(r) is some constant. In the following we denote Cj as some constant that can be 
a function of r. Equation (C.2) follows from the fact that for a 0 ws have that Pija —>■ 
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Cl < 1, z = 1,2, fis+xjoL —>■ C 2 [in Lemma 5.15(i)(a) and (i)(c)]; h{r) j ^ —>■ C 3 (r) 

[m Appendix B, part (c)]; i_{a, fis+i,r) —>■ C 4 (r) [m Lemma 5.15(i)(d)]; /3s+i —)• 0 [T^ 
Corollary 5.14]; and C 5 (a 0 in (5.46) of (Hj)- 

Finally, for r = 0,1,..., s — 2, 


, A(a, r + 1 ) , A(a, r + l)/a*^ 

lim — -r— = lim-^-;-——^- 

ato A{a,r) 04.0 av® A{a,r) 


lima^^^ 

atO 


c(r + 1 ) 
c(r) 


= 0 


and for r = s — 1 , 


^4(0^ = lim ^_ 

ato A(a,s-1} ato a(®-i)/^ A(a,r)/a(®-i)/^ c(s - 1) 


This concludes the proof. 


(C.3) 


(C.4) 
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